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Numerical schemes for continuum models of reaction-diffusion systems 

subject to internal noise 
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We present new numerical schemes to integrate stochastic partial differential equations which 
describe the spatio-temporal dynamics of reaction-diffusion (RD) problems under the effect of in- 
ternal fluctuations. The schemes conserve the nonnegativity of the solutions and incorporate the 
Poissonian nature of internal fluctuations at small densities, their performance being limited by the 
level of approximation of density fluctuations at small scales. We apply the new schemes to two 
different aspects of the Reggeon model namely, the study of its non-equilibrium phase transition 
and the dynamics of fluctuating pulled fronts. In the latter case, our approach allows to reproduce 
quantitatively for the first time microscopic properties within the continuum model. 

PACS numbers: 05.40.-a,05.70.Ln,68.35.Ct 



Continuum representations of the dynamics of 
spatially-extended systems subject to fluctuations is a 
very active area of research in statistical mechanics and 
nonlinear dynamics 0, 0, 0, 0, @ ■ This is because they 
are frequently more tractable than discrete models, they 
can be put forward using simple symmetry arguments 
and applying conservation laws, and therefore they pro- 
vide minimal representations of the observed phenom- 
ena. Important instances are Langevin equations for the 
relaxational dynamics of equilibrium models y|, growth 
interface phenomena Q or coarse-grained descriptions of 
microscopic RD problems 0, |(| . Despite their apparent 
simplicity, most of these models can not be solved ana- 
lytically and one has to resort to approximate analytical 
techniques, or to numerical integration of the stochastic 
time-dependent set of equations using well established 
algorithms Q. In the important instance of RD sys- 
tems subject to internal fluctuations the configurations 
are given by a non-negative density field p{x, t) subject to 
fluctuations of typical strength yj p{x, t) which accounts 
for the Poissonian fluctuations of the number of parti- 
cles at x Unfortunately, standard algorithms fail to 
guarantee both the essential non-negativity of p(x, t) and 
the Poissonian character of its fluctuations. Our purpose 
in this paper is to propose efficient numerical algorithms 
to overcome these problems which will allow us to prove 
the importance of internal fluctuations and to check the 
relevance of their correct description at different scales. 

In this paper we concentrate in the so called Reggeon 
model, which in one dimension is given by Q 



757 = D fc2 +P~P 2 + V^PV(x,t), (1) 



dp 
dt 



where r)(x,t) is a Gaussian white noise. The Reggeon 
model can be obtained under some approximations from 
the microscopic Master equations of RD microscopic 
models using well-known techniques [fj 0]. Heuristi- 
cally, Eq. Q can also be considered as the simplest dy- 
namical equation for a coarse-grained density field with 



a = 1/N, N being the mean-field number of particles 
per site. The Reggeon model provides a minimal repre- 
sentation of the Directed Percolation (DP) universality 
class, which is currently regarded as paradigm of non- 
equilibrium systems with absorbing states jfj: if p(t) is 
the mean density spatial average, there exists a critical 
value of a for which Q undergoes a transition between 
an active phase lim^oo p(t) ^ and an absorbing phase 
for which limt_ >00 p(t) = 0. 

In addition, when a — Eq. (JIJ becomes the so called 
Fishcr-Kolmogorov-Petrovsky-Piscounov (FKPP) equa- 
tion |9( , which displays pulled fronts in which the active 
phase invades the absorbing state flfllllllTp. . Simu- 
lations of microscopic particle models [lllllj] have shown 
that the dynamics of pulled fronts are extremely sensitive 
to microscopic fluctuations at p ~ 1/N, leading to strong 
corrections in the front properties when compared with 
those of the FKPP equation. Since Eq. is usually held 
as a continuum description of some particle models at the 
mesoscopic level (i.e. when p ^> 1/N) one might doubt 
that the Reggeon model describes correctly the behavior 
of pulled fronts subject to internal fluctuations. The ef- 
ficiency and accuracy of the numerical schemes proposed 
here will allow us to show that Eq. indeed incorpo- 
rates the ingredients to explain (even quantitatively) the 
phenomena observed in particle models, thus providing 
also a minimal representation of pulled fronts subject to 
internal fluctuations. 

To simplify the discussion, let us consider the simplest 
possible case for the dynamics of a density subject to 
internal fluctuations: 



dp 
It 



— = ap + y/apr){t). 



(2) 



Typical explicit or implicit methods based on stochas- 
tic Taylor approximations of @ immediately run into 
problems, since they do not conserve the nonnegativity 
of p(t). For example, the Euler approximation is [7| 



Pt+At = pt + apt At 



(3) 



2 



where AWt are random Gaussian numbers with zero 
mean and At variance. Thus, there is a finite probability 
that pt+At becomes negative, and the numerical integra- 
tion comes to a halt. In order to overcome this problem, 
Dickman proposed an interesting solution based on the 
Euler scheme © and the discretization of the possible 
values of pt as multiples of p* = 0(aAt) [l^. Despite its 
success in reproducing the universality class exponents 
of DP using Q and its application to other situations 
Dickman's algorithm is not really a numerical in- 
tegration of a continuum model. Moreover, no general 
study of its convergence and applicability for other situ- 
ations has been done yet. A more technical solution was 
proposed by Schurz and coworkers using Balanced 
Implicit Methods (BIM), in which implicit Euler meth- 
ods are used to impose the nonnegativity of the solution. 
In the case of Eq. (0) the BIM scheme reads [lif 



Pt+At 



Pt + Atpt + y^pl (AW t + \AW t \) 
l + ^Jp~t\AW t \ 



(4) 



which explicitly implements the constraint pt+At > 0, 
and reduces to the Euler algorithm © up to order O(At) 
[r?l . The BIM scheme is known to have the same or- 
der of convergence as the Euler algorithm, namely, the 
error is 0(s/ At) for approximations of individual trajec- 
tories and O(At) for moments of p(t) 0E3- 

Another approach was taken by Pechenik and Levine 
|l2l | employing the exact conditional probability density 
(CDF) V(pt+At\pt) for the stochastic process satisfy- 
ing J5J, which has been known for some time in econ- 
omy as the Cox-Ingersoll-Ross process [2(J- The CDF 
can be expressed in terms of modified Bessel functions 
and, although it can be sampled numerically using re- 
jection or transformation methods \l2\ . it is computa- 
tionally expensive. Here we propose a more efficient 
procedure, which is based on the following: if we de- 
fine rd{t) = X)iLi x ?W' wne re Xi(t) satisfies dxi/dt = 
axi/2+(a/4) 1 ^ 2 r]i{t) with r]i(t) independent white noises, 
then drd/dt = da/4 + ard + (ard) 1 ^ 2 ^^) which coin- 
cides with (J2J in the limit d — ► 0. Since the equation 
for each Xi(t) is linear, rd(t) is the sum of squares of 
Gaussian random numbers with non-zero mean. Thus 
its probability distribution is related to the x 2 distri- 
bution with d degrees of freedom [2lJ. Specifically, we 
find that p t+ At = r (t + At) = x'oW/( 2k ) where 



2a/[a(e 



aAt 



1)], A = ke a ^ P t and X 'oW 



dom number with a noncentral x 2 distribution with zero 
degrees of freedom and noncentrality parameter A whose 
cumulative distribution function is given by [TR HR Elf 



v[xiw<x]=j: { - x/2) 

3=1 



3 e -A/2 



J- 



Hx 2 



2j 



< x]+e- x ' 2 Q{x), 



where x% is a X 2 random number with 2j degrees of free- 
dom and Q(x) is the step function. Equation JSJ is im- 
portant for two reasons: (i) it shows that there is a finite 
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FIG. 1: Conditional probability density as a function of pt+At 
for Eq. J2J with At = 0.1, p t = 1 (left) and p t = 1CT 2 (right). 
Solid lines are the exact solution from JSJ, while dashed lines 
are the approximations obtained using the Euler scheme JHJ. 



probability V(pt+At = 0) = e~ A / 2 for getting into the ab- 
sorbing state, and more importantly (ii) it reveals that 
the probability distribution of x'o W * s a linear combina- 
tion of x 2 probability distributions with Poisson weights. 
This fact can be exploited to generate Pt+At efficiently: 
if we choose K from a Poisson distribution with mean 
A/2, then 



Pt+At 



1 

2k 



if K = 0, 

■ 1 ■ 



if K f 0, 



(0) 



where Zi are independent Gaussian random numbers with 
zero mean and unit variance. 

Another interesting feature of the exact CDF for @ is 
the fact that it converges asymptotically towards the Eu- 
ler approximation (0) when A ~ p t / (aAt) 1 0, . 
However, for small A, the Euler approximation under- 
estimates the large fluctuations present in the exact so- 
lution of J2J). This effect, which can be seen in Fig. ^ 
is related to the fact that the Gaussian approximation 
© of a Poisson random number is only valid when 
the mean value is large enough [19f . The failure of ap- 
proximations like J21 or Q to reproduce large density 
fluctuations at small values of A introduces an effective 
microscopic cutoff p* = 0(aAt) in the numerical simu- 
lations below which these approximations break down. 

Although the scheme (@J can be easily generalized to 
integrate equations like Q), this is not the case for the 
exact sampling of the CDF for J2J. Thus, a splittinq- 
step strategy for integrating Eq. (Q) was proposed in [12| , 
where the time interval At is split into two steps: (i) given 
p t , we use © to integrate (0) and get an intermediate 
value pt+At! (ii) we take pt+At as the initial condition for 
(5) dp/dt = d 2 p/dx 2 — p 2 , producing pt+At with the aid of 



any deterministic numerical algorithm. It can be proved 
that this splitting step method (SSM) converges towards 
the solutions of its order of convergence being O(At) 
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FIG. 2: Left: convergence analysis for the different algo- 
rithms. Points are the critical value of a as a function At 
(below) and of At 1 / 2 (up) while dashed lines are linear fits to 
the data. System size is L = 400. Right: time dependence 
of the mean density p(t) at the critical point a = a c (At) 
with At = 10 -2 obtained using the different algorithms (lines 
are shifted vertically for clarity). Thin line is the power law 
p(t) ~ t' s with S = 0.1595. System size is L = 1000. 



both for realizations and for moments of p{t) [l9(. This 
means that the splitting-step method provides better ap- 
proximations than those based on Euler methods (like 
the Dickman and BIM algorithms) for any realization of 
the noise. This has significant consequences when charac- 
terizing the critical point, as will be shown below. In the 
following we apply the two methods proposed here [BIM 
and the SSM using JfJJl] and the Dickman algorithm to 
the two problems for which Eq. JIJ is archetypal 22]. 

Study of the DP phase transition. To test the proposed 
algorithms, we study the well known non-equilibrium 
phase transition that Reggeon model displays for mod- 
erate values of a 0|. At the critical point, the mean 
average density p(t) = i- J2 X (p{ x > f )) decays like a power 



law p ~ t~ d with 6 ~ 0.1595 |6|. As in [lfj, we identify 
the critical point as the value of a for which we observe 
such a power law decay in p(t). Results for the different 
algorithms are shown in Fig. El where we report the value 
of o~ c as a function of the time step At. As expected, the 
order of convergence of the Dickman and BIM methods is 
0(\/AT), while the SSM has O(At) order of convergence. 
The improvement in the order of convergence comes with 
a price: the computer time needed for our numerical sim- 
ulations at the critical point (see table indicate that 
methods based on Euler approximations, despite having 
an effective microscopic cutoff at p* — 0(oAt), are faster 
than the SSM, and thus could provide better strategies 
for integrating numerically equations for RD models close 
to the critical point, where only accurate approximations 
of large length and time scales are needed. 

Dynamics of fluctuating pulled fronts. When a = 0, 
equation Q displays a wave-like solution (front) which 
travels with velocity vq = 2\fD (provided sharp enough 



TABLE I: Comparison of CPU time at a = <r c (At) with 
At = 10" 2 and L = 400 for the different algorithms in Fig.0 
normalized to that of Dickman's algorithm. 



Method 




Trun 


Dickman 


2.55 


l 


BIM 


1.61 


1.2 


Splitting-Step 


1.45 


7.6 



initial conditions are given) 0, ^} . The dynamics of this 
pulled front is severely affected when microscopic fluctu- 
ations close to the absorbing state p = are considered. 
Specifically, it has been observed in particle models whose 
mean field limit is given by the FKPP equation, that the 
front speed is universally modified as |10L IllL Il4j 



vn = lim 

t — >oc 



•(*)> 



c 



t 



vo - vo 



In 2 AT' 



(7) 



where Xf(t) is the instantaneous position of the front, C 
is a positive constant and N is the number of particles 
per site 0|- Moreover, the pulled front diffuses with 
diffusion constant 



D 



f,N 



= lim 

t — >oc 



((x f (t)-v N t) 2 ) 



C 



2t 



ln 3 A' 



(8) 



where C is a positive constant. Whereas the velocity 
correction can be easily understood because microscopic 
fluctuations at p ~ A -1 provide an effective cutoff in the 
dynamics 0, the diffusion coefficient seems to depend 
on the existence of relatively large fluctuations in the 
density at p ~ A -1 and on their slow relaxation by the 
pulled front dynamics [Til ]. 

As mentioned in the introduction, one might doubt 
that the large microscopic fluctuations at p = A -1 ob- 
served in particle models are correctly reproduced by 
such type of equation like Q. Note, however, that the 
relationship between particle models and the Reggeon 
field model is deep er than at the coarse-grained level. 
Specifically, in 13] it was shown that there is an exact 
duality transformation between the A <-> A + A micro- 
scopic particle model and the so called stochastic FKPP 
equation, which is similar to the Reggeon model but 
with a \/crp(l — p) J](x,t) noise term. For d < 1, the 
noise is only relevant at very small values of p where 
\Jap(\ — p) ~ \fo~p and thus, both the Reggeon model 
and the stochastic FKPP should provide similar results. 

Our results for the front diffusion coefficient, obtained 
by numerical integration of Eq. are reported in Fig. 
Altogether with those of hybrid Monte-Carlo results for 
the A <->■ A + A particle model • ^ s we can see i f° r a 
given time step At, the SSM reproduces the ln~ 3 A re- 
sults for particle models © thus confirming the duality 
relationship between the particle model and the contin- 
uum equation even at the quantitative level. However, 
the other algorithms are more consistent with a ln~ 6 A 



4 




FIG. 3: Front diffusion coefficient as a function of a for the 
different algorithms and different At, compared with hybrid 
MC simulations of the microscopic model A <-> A + A p^|. 
Solid (dashed) line is the log" 3 N (log -6 N) power law. 

scaling which, interestingly, can be obtained through 
standard perturbation techniques based on Gaussian ap- 
proximations for the fluctuations of the front position 
|l0| . The reason for this difference among the various 
schemes is related to the fact that both the Dickman and 
BIM algorithms are based on Gaussian approximations 
for the density fluctuations which are underestimated for 
p < p* = 0(aAt), while only the SSM reproduces exactly 
the large density Poissonian fluctuations (also observed 
in particle models) when p is small. This does not mean 
that the Dickman and BIM algorithms do not converge 
in this case: specifically, if we take At — > we observe 
that the value of the diffusion coefficient approaches that 
of the hybrid MC simulations for the A <-> A+ A (see Fig. 
OJ) . Thus the applicability of the Dickman and BIM algo- 
rithms is limited in this case since they fail to reproduce 
fluctuations at small density and time scales. 

In summary, we have presented new strategies for in- 
tegrating stochastic (partial) differential equations for 
models of RD subject to internal fluctuations. While 
all of them preserve the nonnegativity of the solution, 
algorithms based on Gaussian approximations introduce 
a microscopic cutoff below which density fluctuations are 
not correctly accounted for. This is not important when 
the system properties are dominated by the dynamics of 
large length and time scales (as in critical behavior) , and 
thus, schemes based on Euler approximation suffice to 
integrate numerically equations like l[T]l. However, when 
the observed phenomena are sensitive to microscopic fluc- 
tuations, only algorithms which take into account the 
exact sampling of density fluctuations at small scales 
are computationally efficient. Moreover, our results vali- 
date continuum models like JQ) to study the dynamics of 
fluctuating pulled fronts and corroborate the importance 
of Poissonian large fluctuations of the density at small 
scales. We hope that our results will be used in future 



for the analytical understanding of pulled front dynamics 

Finally, we mention that the methods presented here 
(the BIM and SSM) can be easily extended to other 
situations in which the relevant degrees of freedom are 
non-negative |15i llflj . like the study of density fluctu- 
tions in more general RD problems the understand- 
ing of critical phenomena of systems subject to exter- 
nal/multiplicative noise (e.g. with pr](x,t) noises) 
or the nonlinear modelling of the behavior of interest 
rates in economy ^1 ■ 
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